################################################################################
# GLMMARCodeVersion2 ---this is a revised version of GLMMARCode               
#                  The major modification is to add the option of which        
#                  chain is to be monitored in the MCMC simulation             
#                  Another big change is that this can support the models      
#                  with only random-intercepts---without group-level predictor 
#                  or without random-specific covariates
#                  it has a back-up version : GLMMARCode(revised)       
#                                                                              
# Author: Xun Pang                                                             
#                                                                              
# Modified Date: 1/1/2009                                                      
#                1/2/2009                                                      
#                                                                              
# Usage:  This program supports the following models                           
#          1, binary/ordinal GLMM with two varing-intercepts: time and unit    
#          2, binary/ordinal GLMM with unit/time-specific effects but without  
#               group-level predictors                                         
#          3, binary/ordinal GLMM with two varing-intercepts, but one/both     
#               intercepts are explain by group-level predictors               
#          4, binary/ordinal GLMM with random effects (covariates more than 1  
#              and group-level predictors                                      
#          5, all above models with/without serial correlation                 
#                                                                              
# Program Test Summary:                                                        
#   1/2/2009: the program has been briefly tested for the following            
#             model specifications:                                            
#             random intercepts (two) without group level predictors           
#             random intercpets (two) with one group having predictors         
#             random intercepts (two) with two groups of predictors            
#             one random intercepts, one group of specific random effects      
#                    with/without one/two groups of predictor                  
#            two groups of specific random effects                             
#               with/without one/two groups of pre                             
#             At this time, the program does not support one clustering
#
#  Suggestions for the future: the data arrangement function should be taken out
#           and support the specification that different group-level predictors
#           for different random effects
################################################################################
              
############################ Code Begins #######################################

#________________ *I: Loading necessary packages and set enviroment *__________#

library(MCMCpack);  library(mvtnorm);  library(lattice); library(msm); library(MASS); library(agce); library(Matrix); library(corpcor); library(accuracy); library(bayesSurv); library(smoothSurv); library(bayesm)

#____________________*  Part II: Functions for MCMC Simulation *_____________________#
######################################################################################
#   Instructions for auguments:                                                      
#    y              response variables: binary or ordinal                            
#    Fixed          covariates having unmodeled effects                              
#    UnitRandom     covariates having modeled effects (unit grouping)                
#                   if only a constant, assigne a vector of 1                        
#    TimeRandom     covariates having modeled effects (time grouping)                
#                   if only a constant, assigne a vector of 1                        
#    UnitPred       unit-level predictors                                            
#                   if no group-level predictor, assigne"NULL"                       
#    TimePred       time-level predictors,                                           
#                   if no group-level predictor, assigne"NULL"                       
#    Unit           unit-indicators: numerical                                       
#    Time           time-indicators: numerical                                       
#    m              MCMC iterations                                                  
#    burnin         burnin                                                           
#    init...        initial values                                                   
#    ...0           priors                                                           
#                   whenever W or S includes two or more covariates, D0 and E0       
#                     are matrices                                                   
#    lag            number of lags for autoregressive errors: value permitted: 1, 2  
#    pgm            1-- using PGM-MGMC updating; 0--don't use                        
#    thin           thin the chain of ystar and u                                    
#    monitor        parameters to be monitored and returned. If W or S only include  
#                   random intercetps, betai/betat should not be monitored           
######################################################################################

#////////////////////// CODE 1: Binary or Ordinal without Serial Correlation ////////#

# Starting Values: defaults are the MLE which will be provided in the job files

CholPGM.BTT <- function (y, Fixed, UnitRandom, TimeRandom, UnitPred="NULL", TimePred="NULL", Unit, Time, m=10000, burnin=5000, a1=3, a2=0.1,  beta0=rep(0,24),B0=diag(10, 24), D0=diag(1.67, 5) , d0=9, E0=diag(1.67, 2), e0=6, init.tau, cut0.def, init.ystar=0, init.beta=0,init.b=0, init.c=0,  pgm=1, unitint.add=0, timeint.add=0, monitor=c("rho", "beta", "bi", "ct", "D", "E", "ystar",  "tau", "betai", "betat"), thin=10){

#***************************** Arrange the Data *************************************#
# With/without group predictors results in different data arrangement
# Also, if group predictor only includes a intercept, this intercept will be dropped if 
# at the higher level, there is only a random intercept included
  NewData <- RegressionData(yob=y, Unit=Unit, Time=Time, Fixed=Fixed, UnitRandom=UnitRandom, TimeRandom=TimeRandom, UnitPred=UnitPred, TimePred=TimePred)

#Get results from the approprate function
 NewData <- RegressionData(yob=y, Unit=Unit, Time=Time, Fixed=Fixed, UnitRandom=UnitRandom, TimeRandom=TimeRandom, UnitPred=UnitPred, TimePred=TimePred)

#Get results from the approprate function
  newid1 <- NewData[[1]]
  newid2 <- NewData[[2]]
  Y <- NewData[[3]]
  X <- as.matrix(NewData[[4]])
  W <- as.matrix(NewData[[5]])
  if (unitint.add==1){
    W <- cbind(1, W)
  }
  S <- as.matrix(NewData[[6]])
  if (timeint.add==1){
    S <- cbind(1, S)
  } 
#****************************** Dimensions  ************************************#
   
   # number of all observations
   N <- length(Y)
   # Dimension of all fixed parameter vector beta 
   FB <- ncol(X)
   # Dimension of random fixed parameter vector bi
   RB <- ncol(W)
   # Dimension of random fixed parameter vector ci
   TB <- ncol(S)
 
   # create indices for two non-nested groupings ##########
     entry <- seq(1: N)                   # the global length of observations
     id1 <- Index(newid1)                 # unit id repeating t times t=# observations
     id2 <- (newid2)-min(newid2)+1        # year id repeating n times n=# unit in each time period
     GU <- max(id1)                       # how many units
     GT <- max(id2)                       # how many time periods in total
     categ <- length(unique(Y))           # binary or ordinal
     cat("binary? ",  categ,"\n", sep=" ")
    
   # output the summary of observations
      sample.size.unit <- as.vector (table (id1))
      sample.size.time <- as.vector (table (id2))
      cat("number of observations for each unit: ",  sample.size.unit,"\n", spe=" ")
      cat("number of observations for each time period: ",  sample.size.time,"\n", sep=" ")
    
#_______________________initial values ________________________________________#

      beta <- matrix(init.beta, nrow=1, ncol=FB)
  # This is a big matrix (vectorized), there are many parameters: each unit has RB random effects
      bbb <- matrix(init.b, nrow=1, ncol=RB*GU)
  # The Wishart random generator produces the lower triangular of the draws
      DDD <- matrix(0, nrow=1, ncol=(RB^2-RB)/2+RB)
  # matrix with colomn = #time periods * time periods
      ccc <- matrix(init.c, nrow=1, ncol=GT*TB)
# The Wishart random generator produces the lower triangular of the draws
      EEE <- matrix(0, nrow=1, ncol=(TB^2-TB)/2+TB) 
      betaii <- matrix(0, nrow=1, ncol=RB*GU)
      gammai <- matrix (0, nrow=1, ncol=TB*GT)
      ystarm <- matrix(init.ystar, nrow=1, ncol=N)
        if (categ>2){
        tau <- matrix(init.tau, nrow=1, ncol=(categ-2))
          }
   
 #**************************** loop begins *************************************#
 #Tracking the process of MCMC simulation
   
    for(g in 2:(m+burnin)) { 
       if ((g%%5 ==0)| g==m){
          cat("g=",g,"\n")
        }
       
#________________ Values of the previous iteration _____________________________
# use ~[1,] instead of init.~ because of the dimension confusion
       
   if (g==2){
   current.bi <- matrix(bbb[1,], ncol=RB)           # current values of bi
   current.ct <- matrix(ccc[1,], ncol=TB)           # current values of ct
   current.beta <- beta[1,]                         # current values of beta
   current.ystar <- ystarm[1,]
   if (categ>2){
     current.tau <- tau[1,]
   }
 }else{
   current.bi <- new.bi                            # current values of bi
   current.ct <- new.ct                            # current values of ct
   current.beta <- new.beta                        # current values of beta
   current.ystar <- new.ystar                      # the current value of ystar
  if (categ>2){
     current.tau <- new.tau
   }
 }

#***************************  Block 1: variance Matrix of E *******************#
       
       if (ncol(S)==1){                       # if S only contains 1 element (most likely a constant)
         Ei<-rgamma(1,(e0+GT)/2, rate=(E0+sum(current.ct^2))/2)
         Ei <- as.matrix(Ei)
          }else{
        OEin <- solve(E0)+sumtrans(Mat=current.ct)
        sE<-rWishart (1, e0+GT, solve(OEin))
        Ei <- Wishart(element=sE, col=TB)
        }
# Block 1 end
       
 #***************************  Block 2: variance Matrix of D *******************#

    if (ncol(W)==1){                         # if W only contains 1 element (most likely a constant)
         Di<-rgamma(1,(d0+GT)/2, rate=(D0+sum(current.bi^2))/2)
         Di <- as.matrix(Di)
          }else{
       ODin <- solve(D0)+sumtrans(Mat=current.bi)
       sD<-rWishart (1, d0+GU, solve(ODin))
       Di <- Wishart(element=sD, col=RB)
         }
# Block 2 end

#####################################################################################
# Define variables for outputs from the loop 

      betaterm1 <- matrix(0, nrow=FB, ncol=FB)    # for posterior covariance of beta
      betaterm2 <- rep(0, FB)                     # for posterior mean of beta
      new.ystar <- rep(NA, N)                    # holder for next iteration of ystar

# Big loop for data augmentation and values for beta 

   for (n in 1: GU){                              
      rows <- entry[which(id1==n)]                    # the obserbations of unit n
      nT <- length(rows)                              # length of those observations
      yi <- Y[rows]                                   # Observed responses of n
      ystari <- current.ystar[rows]                   # current latent responses of n
      Xi <- as.matrix(X[rows,])
      Wi <- as.matrix(W[rows,])
      Si <- as.matrix(S[rows,])
      bbbn <- current.bi[n,]
      TT <- id2[rows]
      Tccin <- as.matrix(current.ct[TT,])

 # the term of Si*c_{T_i}
       rant <- c()
      for (l in 1:nT){                                # Using loop to do it
          rant[l] <- Si[l,]%*%cbind(Tccin[l,])
            }
      
#******************************* Block 3: Data Augmentation ********************#
        mean.ystar <- Xi%*%cbind(current.beta) +Wi%*% cbind(bbbn)+rant
        ystar <- c()
     if (categ==2){
        for (r in 1:nT){
         ystar[r] <- BTnormlog(x=yi[r], mu=mean.ystar[r], sigma=1)
           }
              }else{
         for (r in 1:nT){
         ystar[r] <- OrdinalAugment(obs=yi[r], mu=mean.ystar[r], sigma=1, maxcat=max(Y), mincat=min(Y), cut0=cut0.def, cutpoint=current.tau)
             }
      }
  
        new.ystar[rows] <- ystar
     

 # calculate BB1 and BB2 for beta 
    
       var.beta <- pseudoinverse(diag(nT)+Wi%*%solve(Di)%*%t(Wi))
       betaterm1 <- betaterm1+t(Xi) %*% var.beta %*%Xi
       betaterm2 <- betaterm2+t(Xi)%*% var.beta %*% (ystar-rant)
         }
      
 #***************************Block 4 Beta and {b_i} ************************************#
       # 1) updating beta
       B1<-solve(betaterm1+solve(B0))
       beta1<-t(B1%*%(betaterm2+solve(B0)%*%beta0))
      # B1 <- make.positive.definite(B1)
       new.beta<-mvrnorm(1, beta1, B1)

      # 2)The updating  {b_i}
       
      new.bi <- matrix(NA, nrow=1, ncol=RB)

    for (n in 1: GU){                                   # loop for each unit
        rows <- entry[which(id1==n)]                    # the obserbations of unit n
        nT <- length(rows)
        ystari <- new.ystar[rows]
        Xi <- as.matrix(X[rows,])
        Wi <- as.matrix(W[rows,])
        Si <- as.matrix(S[rows,])
        TT <- id2[rows]
        Tccin <- as.matrix(current.ct[TT,])
     
         rant <- c()
         for (l in 1:nT){
          rant[l] <- Si[l,]%*%cbind(Tccin[l,])
          }
     
      Cov.b <- solve(Di + t(Wi) %*%Wi)
      mean.b <- Cov.b%*% t(Wi) %*% (ystari-Xi%*% cbind(new.beta)-rant)
      if (length(mean.b)==1){
         b.drawn <- as.matrix(rnorm(1, mean.b, sqrt(Cov.b)))
       }else{
        # Cov.b <- make.positive.definite(Cov.b)
         b.drawn <- mvrnorm(1, mu=mean.b, Sigma=Cov.b)
         }
          new.bi <- rbind(new.bi, t(b.drawn))
        }
       
          new.bi <- as.matrix(new.bi[-1,])
           
#************************** Block 5: Time-Specific Random Effect ****************#

         new.ct <-matrix(NA, nrow=1, ncol=TB)
         for (j in 1: GT){
           row1 <- entry[which(id2==j)]    # Those entries with t=j
           H <- length(row1)               #  #observations at t period
           iden1 <- id1[row1]              # The units at t period
           tttn <- as.matrix(new.bi[iden1,])         # The random unit-specific effects of
           yt1 <- new.ystar[row1]
           Xt1 <- as.matrix(X[row1,])
           Wt1<- as.matrix(W[row1,])
           St1 <- as.matrix(S[row1,])
            
            term7 <- c()
             for (q in 1:H){
              term7[q] <- Wt1[q,]%*% cbind(tttn[q,])
              }

           cov.S <- solve(Ei+t(St1)%*%St1)
           mean.S <- cov.S%*% (t(St1) %*%(yt1-term7-Xt1%*%cbind(new.beta)))
             if (length(mean.S)==1){
               supd <- as.matrix(rnorm(1, mean.S, sqrt(cov.S)))
               } else {
               supd <- rmvnorm(1, mean.S, cov.S)
               }
                 new.ct <- rbind(new.ct, supd)
           }
       
            new.ct <- as.matrix(new.ct[-1,])

#******************************Block6:  Cutpoints Updating ******************#

       if (categ>2){
       new.tau <- CutPointUpdate(obs=Y, latent=new.ystar, cut0=cut0.def, oldcut=current.tau)
          }

# **************************** PGM Updating *******************************#
       
       if (pgm==1){                    # Make choice whether to use updating (1 yes)
           if (any(beta0!=0)){
             warning("priors have to be centered at 0 for using PGM-MGMC updating")
             return(NULL)
             }
        
            qq1 <- 0
            qq2 <- 0
          for (f in 1: GU){
            row2 <- entry[which(id1==f)]    # Those entries with n=i
            H <- length(row2)               
            iden3 <-id2[row2]               # The units at t period
            tttn <- new.bi[f,]                 # The random unit-specific effects
            yt1 <- new.ystar[row2]
            Xt1 <- as.matrix(X[row2,])
            Wt1 <- as.matrix(W[row2,])
            St1 <- as.matrix(S[row2,])
            
            timec <- as.matrix(new.ct[iden3,])
               rrant <- c()
              for (w in 1:H){
                rrant[w] <- St1[w,]%*% cbind(timec[w,])
                }
            
             termm <- yt1-Xt1%*%as.matrix(new.beta)-Wt1%*%as.matrix(tttn)-rrant
             qq1 <- qq1+t(termm)%*%termm
             qq2 <- qq2+rbind(tttn)%*%Di%*%cbind(tttn)
            }
           
             qq3 <- rbind(new.beta)%*% solve(B0)%*% cbind(new.beta)
           
             qq4 <- 0
             for (t in 1: GT){
               qq4 <- qq4+rbind(new.ct[t,])%*%Ei%*%cbind(new.ct[t,])
               }
           
             bb <- (qq1+qq2+qq3+qq4)/2
                         
           if (categ>2){
               aa <- (N+ ncol(beta)+ncol(bbb)+ncol(ccc)+length(new.tau)+1)/2
                  } else {
               aa <- (N+ ncol(beta)+ncol(bbb)+ncol(ccc)+1)/2
           }

           gamm <- sqrt(rgamma(1, shape=aa, scale=1/bb))
           new.ystar <- new.ystar*gamm
           new.beta <- new.beta*gamm
           new.bi <- new.bi*gamm
           new.ct <- new.ct*gamm
           
           if (categ>2){
              new.tau <- new.tau*gamm
           }
             cat("gamm ", gamm, "\n", sep=" ")
         }
       
#*********************** Updata ******************************************
# Compute time- and unit-specific random effects
# if only random intercept is specified, there is no specific random effects
#    will be computed
# Be very careful about the dimensions
# 
 if ("betai"%in% monitor){
    if (length(unique(W))>1){                     # whether only random intercept?
       random.pred <- as.matrix(new.bi)           # the new draw of random unit effects

     # get the group-level average effects
    if (length(unique(S))==1){                    # whether S is only constant?
      if (!is.numeric(TimePred)){                 # if S is constant and no group-level predictor
       beta.pp <- new.beta[(ncol(Fixed)+1):(ncol(X))]  # Average effects at the unit-level
        }else{          #if S is constant, but group-level predictors
       beta.pp <- new.beta[(ncol(Fixed)+1):(ncol(X)-ncol(TimePred))] 
     }
    }
    if (length(unique(S))>1){              # if S contains covariates other than a constant
         if(is.numeric(TimePred)){         # if S and group level predictors
        beta.pp <- new.beta[(ncol(Fixed)+1):(ncol(X)-ncol(TimePred)*ncol(S))]
      } else{                              # if no group level predictors
         beta.pp <- new.beta[(ncol(Fixed)+1):(ncol(X)-ncol(S))]
       }
       }
       
     # Compute unit-specific random effects    
    if (is.numeric(UnitPred)){             # if group level predictors at the unit level
       beta.pred <- matrix(beta.pp, ncol=ncol(UnitPred), byrow=TRUE)  # matrix form of the average effects
       UnitCova <- remove.NA(unique(UnitPred))                        # get rid of missing data
       beta.2 <- Random.Coef(Cova=UnitCova, random=random.pred, coef.pred=beta.pred, RBB=RB )
       betaii <- rbind(betaii, beta.2)
 }else{                                    # if no group level predictors
      beta.2 <- matrix(NA, ncol=ncol(random.pred), nrow=GU)  # simply add the random residual on the average
      for (h in 1: GU){
        beta.2[h,] <- random.pred[h,]+beta.pp
      }
      betaii <- rbind(betaii, as.vector(beta.2))
    } 
  } else {                                 # if only random intercept specified, stop grogram and return error message
       warning("No unit-specific effects to return")
      return (NULL)
     }
  }

# Compute time-specific random effects
  if  ("betat"%in% monitor){
       time.random <- as.matrix(new.ct)
      if (length(unique(S))>1){
        if (is.numeric(TimePred)){
          gamma.pp <- new.beta[(ncol(X)-ncol(TimePred)*ncol(S)+1):ncol(X)]
          time.pred <- matrix(gamma.pp, ncol=ncol(TimePred), byrow=TRUE)
          TimeCC <- as.matrix(na.omit((unique(TimePred))))
          gamma.2 <- Random.Coef(Cova=TimeCC, random=time.random, coef.pred=time.pred, RBB=TB)
          gammai <- rbind(gammai, gamma.2)
       }else{
          gamma.pp <- new.beta[(ncol(X)-ncol(S)+1):ncol(X)]
          gamma.2 <- matrix(NA, ncol=ncol(time.random), nrow=GT)
         for (h in 1: GT){
           gamma.2[h,] <- time.random[h,]+gamma.pp
          }
         gammai <- rbind(gammai, as.vector(gamma.2))
       } 
     }else{
        warning("No time-specific effects to return")
        return (NULL)
      }
   }

#********************** Updata other parameters *********************************#
      if ("ystar"%in% monitor){
          if ((g%%thin ==0)| g==m){
         ystarm <- rbind(ystarm, new.ystar)
          }
        }
       
       if ("bi" %in% monitor){
         bbb <- rbind(bbb, as.vector(new.bi))
       }
       
       if ("ct"%in% monitor){
        timetime <- as.vector(new.ct)
        ccc<-rbind(ccc, timetime)
       }
       
       if ("beta"%in% monitor){
        beta<-rbind(beta,new.beta)
       }
       
       if ("D"%in% monitor){
         if (ncol(W)==1){
           DDD <- rbind(DDD, as.vector(Di))
         }else{
           DDD <- rbind(DDD, sD)
         }
       }
       
        if ("E"%in% monitor){
            if (ncol(S)==1){
           EEE <- rbind(EEE, as.vector(Ei))
         }else{
           EEE <- rbind(EEE, sE)
         }
       }
       
       if (categ>2){
        if ("tau"%in% monitor){
        tau <- rbind(tau, new.tau)
       }
     }  
     }

#************************** Returning MCMC output**********************************
      
  if ("beta"%in% monitor){ 
     posterior.beta <- beta[-c(1:burnin),]
  }else{
     posterior.beta <- beta
   }
  if ("bi"%in% monitor){ 
     posterior.unit.res <- bbb[-c(1:burnin),]
  }else{
     posterior.unit.res <- bbb
  }
 if ("ct"%in% monitor){ 
     posterior.time.res <- ccc[-c(1:burnin),]
  }else{
     posterior.time.res <- ccc
  }   
  if ("E"%in% monitor){ 
     posterior.time.cov <- EEE[-c(1:burnin),]
  }else{
     posterior.time.cov <- EEE
  }    
  if ("D"%in% monitor){ 
     posterior.unit.cov <- DDD[-c(1:burnin),]
  }else{
     posterior.unit.cov <- DDD
  }
  if ("betai"%in% monitor){ 
     posterior.unit <- betaii[-c(1:burnin),]
  }else{
     posterior.unit <- betaii
  }   
  if ("betat"%in% monitor){ 
      posterior.time <- gammai[-c(1:burnin),]
  }else{
      posterior.time <- gammai
  }    
   if ("ystar"%in% monitor){ 
     posterior.data <- ystarm[-c(1:(burnin/thin)),]
  }else{
     posterior.data <- ystarm
  }
  if(categ>2){
   if ("tau"%in% monitor){ 
      posterior.tau <- tau[-c(1:burnin),]
  }else{
      posterior.tau <- tau
  }
  } 
       
     if (categ>3){
         return(list( Beta=posterior.beta,  bi=posterior.unit.res, ci=posterior.time.res, Betai=posterior.unit,  Betat=posterior.time , var.c=posterior.time.cov, Cov.betai=posterior.unit.cov, cutpoint=posterior.tau, latentdata=posterior.data))
       }else{
         return(list( Beta=posterior.beta, bi=posterior.unit.res, ct=posterior.time.res , Betai=posterior.unit,  Betat=posterior.time,   var.c=posterior.time.cov, Cov.betai=posterior.unit.cov, latentdata=posterior.data))
       }
}



###################### Function for BTSCS/OTSCS: With modeling serial correlation #######

 CholARP.BTSCSART <- function(y,Fixed, UnitRandom, TimeRandom, UnitPred="NULL", TimePred="NULL", Unit, Time,  m=10000, burnin=5000,  init.tau, init.phi, init.ystar=0, init.beta=0,  init.b=0, init.c=0,a1=3, a2=0.1, beta0=rep(0,24),B0=diag(10, 24), D0=diag(1.67, 5) , d0=9, E0=diag(1.67, 2), e0=6,  pgm=1, lag=2, chol=1,  monitor=c("rho", "beta", "bi", "ct", "D", "E", "ystar", "u", "tau", "betai", "betat"), unitint.add=0, timeint.add=0,thin=100){

#******************* Arrange Data 
# With/without group predictors results in different data arrangement
# Also, if group predictor only includes a intercept, this intercept will be dropped if 
# at the higher level, there is only a random intercept included
  
  NewData <- RegressionData(yob=y, Unit=Unit, Time=Time, Fixed=Fixed, UnitRandom=UnitRandom, TimeRandom=TimeRandom, UnitPred=UnitPred, TimePred=TimePred)

#Get results from the approprate function
  newid1 <- NewData[[1]]
  newid2 <- NewData[[2]]
  Y <- NewData[[3]]
  X <- as.matrix(NewData[[4]])
  W <- as.matrix(NewData[[5]])
  if (unitint.add==1){
   W <- cbind(1, W)
  }
  S <- as.matrix(NewData[[6]])
  if (timeint.add==1){
   S <- cbind(1, S)
  } 
#****************************** Dimensions  ************************************#
   
   # number of all observations
   N <- length(Y)
   cat("N", N, sep=" ")
   # Dimension of all fixed parameter vector beta 
   FB <- ncol(X)
   # Dimension of random fixed parameter vector bi
   RB <- ncol(W)
   # Dimension of random fixed parameter vector ci
   TB <- ncol(S)
 
   # create indices for two non-nested groupings ##########
     entry <- seq(1: N)                   # the global length of observations
     id1 <- Index(newid1)                 # unit id repeating t times t=# observations
     id2 <- (newid2)-min(newid2)+1        # year id repeating n times n=# unit in each time period
     GU <- max(id1)                       # how many units
     GT <- max(id2)                       # how many time periods in total
     categ <- length(unique(Y))           # binary or ordinal
     cat("binary? ",  categ,"\n", sep=" ")
    
   # output the summary of observations
      sample.size.unit <- as.vector (table (id1))
      sample.size.time <- as.vector (table (id2))
      cat("number of observations for each unit: ",  sample.size.unit,"\n", spe=" ")
      cat("number of observations for each time period: ",  sample.size.time,"\n", sep=" ")
    
#**************************** initial values **************************************#
   
     if (lag==1){
        phi <- init.phi
      }else{
        phi <- matrix(init.phi, nrow=1, ncol=lag)
         }
      beta <- matrix(init.beta, nrow=1, ncol=FB)
      bbb <- matrix(init.b, nrow=1, ncol=RB*GU)
      DDD <- matrix(0, nrow=1, ncol=(RB^2-RB)/2+RB)
      ccc <- matrix(init.c, nrow=1, ncol=GT*TB)
      EEE <- matrix(0, nrow=1, ncol=(TB^2-TB)/2+TB)
      betaii <- matrix(0, nrow=1, ncol=RB*GU)
      gammai <- matrix (0, nrow=1, ncol=GT*TB)
      ystarm <- matrix(init.ystar, nrow=1, ncol=N)
      u.aux <-  matrix(rnorm(N), nrow=1, ncol=N)
       if (categ>2){
        tau <- matrix(init.tau, nrow=1, ncol=(categ-2))
          }
   
############### loop begins ################################################

 #Tracking the process of MCMC simulation
    for(g in 2:(m+burnin)) { 
       if ((g%%5 ==0)| g==m){
          cat("g=",g,"\n")
        }
       
# The current values of those parameters

  if (g==2){
   current.bi <- matrix(bbb[1,], ncol=RB)           # current values of bi
   current.ct <- matrix(ccc[1,], ncol=TB)           # current values of ct
   current.beta <- beta[1,]                         # current values of beta
   current.ystar <- ystarm[1,]                      # the current value of ystar
   current.phi <- init.phi
   if (categ>2){
     current.tau <- init.tau
   }
  }else{
   current.bi <- new.bi                            # current values of bi
   current.ct <- new.ct                            # current values of ct
   current.beta <- new.beta                        # current values of beta
   current.ystar <- new.ystar                      # the current value of ystar
   current.phi <- new.phi
  if (categ>2){
     current.tau <- new.tau
   }
 }
  
#_______________Block 1: variance Matrix of E _____________________________
       
      if (ncol(S)==1){                     # if S only contains 1 element (most likely a constant)
         Ei<-rgamma(1,(e0+GT)/2, rate=(E0+sum(current.ct^2))/2)
         Ei <- as.matrix(Ei)
      }else{
        OEin <- solve(E0)+sumtrans(Mat=current.ct)
        sE<-rWishart (1, e0+GT, solve(OEin))
        Ei <- Wishart(element=sE, col=TB)
      }
 # End of Block 1
       
#_______________Block 2: variance Matrix of D _____________________________
       
      if (ncol(W)==1){                        # if W only contains 1 element (most likely a constant)
         Di<-rgamma(1,(d0+GT)/2, rate=(D0+sum(current.bi^2))/2)
         Di <- as.matrix(Di)
      }else{
         ODin <- solve(D0)+sumtrans(Mat=current.bi)
         sD<-rWishart (1, d0+GU, solve(ODin))
         Di <- Wishart(element=sD, col=RB)
      }
  # End of Block 2
 #_______________Compute the variance-covariance matrix of ystar _____
  
    if (lag==1){
      var.phig <- 1/(1- current.phi^2)                             # the variance computed by using the current values of phi
      TCov <- Cov.AR1(col=GT, rho= current.phi,sigma=var.phig)
   }else{
       CCov <- Cov.ARp(phig=current.phi, col=GT)        
       TCov <- CCov[[2]]
       Sig <- TCov[1:lag, 1:lag]
       SSig<-solve(Sig) 
       var <- Sig[1,1]
      }

# Variables for outputs from the loop 
       
      betaterm1 <- matrix(0, nrow=FB, ncol=FB)     # for posterior covariance of beta
      betaterm2 <- rep(0, FB)                      # for posterior mean of beta
      new.u<- rep(NA, N)
      kap <- rep(NA, GU)                           # Store the kappa for variance
      qu <- rep(NA,N)                              # Store the output of the auxiliary term
      new.ystar <- rep(NA, N)

   for (n in 1: GU){                               # loop for each unit
      rows <- entry[which(id1==n)]                 # the obserbations of unit n
      nT <- length(rows)
      yi <- Y[rows]
      ystari <- current.ystar[rows]
      Xi <- as.matrix(X[rows,])
      Wi <- as.matrix(W[rows,])
      Si <- as.matrix(S[rows,])
      bbbn <- current.bi[n,]
      TT <- id2[rows]
      Tccin <- as.matrix(current.ct[TT,])
     
       rant <- c()
         for (l in 1:nT){
          rant[l] <- Si[l,]%*%cbind(Tccin[l,])
        }

   #__________________________Block 3:Updata U _____________________________
     
        nCov <- as.matrix(TCov[1:nT, 1:nT])
        nVK <- V.matrix(Omega=nCov, col=nT)
        nV <- nVK[[2]]                          
        nK <- nVK[[1]]                  
        Cov.U <- solve(diag(nT) + 1/nK*t(nV)%*%nV)   
        mean.U <- Cov.U %*% t(nV) %*% (ystari-Xi%*%cbind(current.beta)-Wi%*%cbind(bbbn)-rant)/nK
        u.draw <- mvrnorm(1, mu=mean.U, Sigma=Cov.U)
        new.u[rows] <- u.draw
        u.decom <- nV%*%cbind(u.draw)
        qu[rows] <- u.decom
        kap[n] <- nK
      
  # End of Block 3
  #__________________________Block 4: Updata ystar ________________________
       if (chol==1){
         mean.ystar <- Xi%*% cbind(current.beta)+Wi%*%cbind(bbbn)+rant+u.decom
         ystar <- c()
           for (r in 1: nT){
             if (categ==2){
                ystar[r] <- BTnormlog(x=yi[r], mu=mean.ystar[r], sigma=sqrt(nK))
             }else{
                ystar[r] <- OrdinalAugment(obs=yi[r], mu=mean.ystar[r], sigma=sqrt(nK), maxcat=max(Y), mincat=min(Y), cutpoint=current.tau)
             }
          }
       }
      
       if (chol==0){
          mean.ystar <-  Xi%*% cbind(current.beta)+rant
          ystar <- DataAug3(yold=current.ystar, ymean=mean.ystar, yi=yi, Cov=(nCov+Wi%*%solve(Di)%*%t(Wi)),  categ=categ, tau=new.tau, Y=Y)
           }
      
       if (chol==2){
          ystar <- DataAug2(Fix=Xi, Rand=Wi, T=T, Dep=yi, Ldep=ystari, FP=current.beta, RP=bbbn, TRP=rant, unit=n, Lag=lag, phig=phig, Sig=Sig, Entry=entry, Id1=id1, Y=Y, tau=current.tau, categ=categ)
        ystar[which(ystar==Inf)] <- 200
         ystar[which(ystar==-Inf)] <- -200
        }
     new.ystar[rows] <- ystar

   # End of Block 4
        
   # calculate BB1 and BB2 for beta
      
      var.beta <- pseudoinverse(nCov+Wi%*%solve(Di)%*%t(Wi))
      betaterm1 <- betaterm1+t(Xi) %*% var.beta %*%Xi
      betaterm2 <- betaterm2+t(Xi)%*% var.beta %*% (ystar-rant)
    }
       
 #_______________Block 5 Beta and {b_i}: ______________________________________
       
       # 1) \beta
       B1<-pseudoinverse(betaterm1+solve(B0))
       beta1<-t(B1%*%(betaterm2+solve(B0)%*%beta0))
       #B1 <- make.positive.definite(B1)
       new.beta<-mvrnorm(1, beta1, B1)
       
       # 2) {b_i}
       new.bi <- matrix(NA, nrow=1, ncol=RB)
       for (n in 1: GU){                                # loop for each unit
         rows <- entry[which(id1==n)]                    # the obserbations of unit n
         nT <- length(rows)
         nCov <- as.matrix(TCov[1:nT, 1:nT])
         ystari <- new.ystar[rows]
         Xi <- as.matrix(X[rows,])
         Wi <- as.matrix(W[rows,])
         Si <- as.matrix(S[rows,])
         TT <- id2[rows]
         Tccin <- as.matrix(current.ct[TT,])
     
         rant <- c()
         for (l in 1:nT){
          rant[l] <- Si[l,]%*%cbind(Tccin[l,])
          }
     
         Cov.b <- solve(Di + t(Wi)%*%solve(nCov)%*%Wi)
         mean.b <- Cov.b%*% t(Wi)%*%solve(nCov)%*%(ystari-Xi%*% cbind(new.beta)-rant)
         if (length(mean.b)==1){
           b.drawn <- as.matrix(rnorm(1, mean.b, sqrt(Cov.b)))
         }else{
          # Cov.b <- make.positive.definite(Cov.b)
           b.drawn <- mvrnorm(1, mu=mean.b, Sigma=Cov.b)
         }
          new.bi <- rbind(new.bi, t(b.drawn))
       }

        new.bi <- as.matrix(new.bi[-1,])

       # End of Block 5
       
  #___________________ Block 6: Time-Specific Random Effect ________________

        new.ct <-matrix(NA, nrow=1, ncol=TB)
        for (j in 1: GT){
            row1<-entry[which(id2==j)]  # Those entries with t=j
            H <- length(row1)           #  #observations at t period
            iden1<-id1[row1]            # The units at t period
            tttn <- as.matrix(new.bi[iden1,])         # The random unit-specific effects of
            yt1 <- new.ystar[row1]
            Xt1 <- as.matrix(X[row1,])
            Wt1<- as.matrix(W[row1,])
            St1 <- as.matrix(S[row1,])
            Ut1 <- qu[row1]
            kapt1 <- kap[iden1]
            
            term7 <- c()
            for (q in 1:H){
               term7[q] <- Wt1[q,]%*% cbind(tttn[q,])
              }
            
            PM <- matrix(0, nrow=H, ncol=H)
            for (c in 1:H){
               PM[c,c] <- kapt1[c]
            }
            
           cov.S <- solve(Ei+t(St1)%*%solve(PM)%*%St1)
           mean.S <- cov.S%*% t(St1) %*%solve(PM)%*%(yt1-term7-Xt1%*%cbind(new.beta)-Ut1)
            if (length(mean.S)==1){
                supd <- as.matrix(rnorm(1, mean.S, sqrt(cov.S)))
             } else {
                supd <- rmvnorm(1, mean.S, cov.S)
             }
          
               new.ct <- rbind(new.ct, supd)
             }
       
           new.ct <- as.matrix(new.ct[-1,])

       # End of Block 6
       
#__________________Block 7: Autoregressive Coef.s_________________________

     if (lag==1){
       PhiPr <- PhiTerm0(Entry=entry, unit=id1, NN=GU, yystar=new.ystar,newb=new.bi, PX=X, PW=W, PS=S, Ti=id2, CC=new.ct, newbeta=new.beta)
       PhiSig <- PhiPr[1]
       PhiMean <- PhiPr[2]
       hatPhi<-1/PhiSig
       hatphi<-hatPhi*PhiMean
       phi.draw <- rtnorm(1, hatphi, sqrt(hatPhi), lower=-1, upper=1)
         cat("phi.draw ", phi.draw, "\n", sep=" ")
        ####### Accept and Reject ########
       
       alpha<-AcceptRate0(rho=phi.draw, orho=current.phi, unit=id1, precision=1,AY=new.ystar, AX=X, AW=W, AS=S, Abeta=new.beta, Ab=new.bi, AT=id2, Ac=new.ct, Entry=entry)

        u <- runif(1)
        if( u <= alpha ) {
            new.phi<-phi.draw
        }else{
           new.phi<-current.phi
        }
    }
       
      if (lag>1){ 
         PhiPr <- PhiTermARpT(Entry=entry, unit=id1, NN=GU, yystar=new.ystar, newb=new.bi, PX=X, PW=W, PS=S, Ti=id2, CC=new.ct, newbeta=new.beta, lag=lag)
         PhiSig <- PhiPr[,-(lag+1)]
         PhiMean <- PhiPr[,(lag+1)]
         phigg<-PhiPropARp(Ymt=PhiSig, Ymt2=PhiMean,scale=1, Oldphi=phig)
         cat("phi.draw ", phigg, "\n", sep=" ")
         
 ####### Accept and Reject ########
         
        alpha<-AcceptRateARpT(nrho=phigg, orho=current.phi, Osig=Sig , unit=id1, AY=new.ystar, AX=X, AW=W, AS=S, Abeta=new.beta, Ab=new.bi,  AT=id2, Ac=new.ct, Entry=entry, f=g, lag=lag)
     #   cat("alpha ", alpha, "\n", sep=" ")
       u <- runif(1)
        if( u <= alpha ) {
            new.phi <- phigg
        } else {
            new.phi<- current.phi
        }
     }
  #________________Block8 _Cutpoints Updating ________________________#

       if (categ>2){
         new.tau <- CutPointUpdate(obs=Y, latent=new.ystar, oldcut=current.tau)
       }

 # __________________________ PGM Updating_________________________________
       
       if (pgm==1){
            if (any(beta0!=0)){
        warning("priors have to be centered at 0 for using PGM-MGMC updating")
        return(NULL)
         }
          phig <-new.phi   
          if (lag==1){
            Phig <- 1/(1-phig^2)
            TCov <- Cov.AR1(col=GT, rho=phig,sigma=Phig)
          }else{
            CCov <- Cov.ARp(phig=phig, col=GT)        
            TCov <- CCov[[2]]
          }
           
          qq1 <- 0
          qq2 <- 0
          for (f in 1: GU){
            row2 <- entry[which(id1==f)]       # Those entries with n=i
            H <- length(row2)               
            iden3 <-id2[row2]                  # The units at t period
            tttn <- new.bi[f,]                 # The random unit-specific effects
            yt1 <- new.ystar[row2]
            Xt1 <- as.matrix(X[row2,])
            Wt1 <- as.matrix(W[row2,])
            St1 <- as.matrix(S[row2,])
            nCov <- as.matrix(TCov[1:H, 1:H])
            timec <- as.matrix(new.ct[iden3,])
            rrant <- c()
            for (w in 1:H){
                rrant[w] <- St1[w,]%*% cbind(timec[w,])
            }
            termm <- yt1-Xt1%*%as.matrix(new.beta)-Wt1%*%as.matrix(tttn)-rrant
            qq1 <- qq1+t(termm)%*%solve(nCov)%*%termm
            qq2 <- qq2+rbind(tttn)%*%Di%*%cbind(tttn)
          }
           
          qq3 <- rbind(new.beta)%*% solve(B0)%*% cbind(new.beta)
           
          qq4 <- 0
          for (t in 1: GT){
            qq4 <- qq4+rbind(new.ct[t,])%*%Ei%*%cbind(new.ct[t,])
          }
           
          bb <- (qq1+qq2+qq3+qq4)/2
                         
          if (categ>2){
               aa <- (N+ ncol(beta)+ncol(bbb)+ncol(ccc)+length(new.tau)+1)/2
          } else{
               aa <- (N+ ncol(beta)+ncol(bbb)+ncol(ccc)+1)/2
          }

          gamm <- sqrt(rgamma(1, shape=aa, scale=1/bb))
          new.ystar <- new.ystar*gamm
          new.beta <- new.beta*gamm
          new.bi <- new.bi*gamm
          new.ct <- new.ct*gamm
           
          if (categ>2){
            new.tau <- new.tau*gamm
          }
            cat("gamm ", gamm, "\n", sep=" ")
          }
       
#*********************** Update the parameters  ******************************************
# Compute time- and unit-specific random effects
# if only random intercept is specified, there is no specific random effects
#    will be computed
# Be very careful about the dimensions
# 
 if ("betai"%in% monitor){
    if (length(unique(W))>1){                     # whether only random intercept?
       random.pred <- as.matrix(new.bi)           # the new draw of random unit effects

     # get the group-level average effects
    if (length(unique(S))==1){                    # whether S is only constant?
      if (!is.numeric(TimePred)){                 # if S is constant and no group-level predictor
          beta.pp <- new.beta[(ncol(Fixed)+1):(ncol(X))]  # Average effects at the unit-level
       }else{          #if S is constant, but group-level predictors
          beta.pp <- new.beta[(ncol(Fixed)+1):(ncol(X)-ncol(TimePred))] 
       }
    }
    if (length(unique(S))>1){              # if S contains covariates other than a constant
       if(is.numeric(TimePred)){         # if S and group level predictors
          beta.pp <- new.beta[(ncol(Fixed)+1):(ncol(X)-ncol(TimePred)*ncol(S))]
       } else{                              # if no group level predictors
          beta.pp <- new.beta[(ncol(Fixed)+1):(ncol(X)-ncol(S))]
       }
    }
       
     # Compute unit-specific random effects    
    if (is.numeric(UnitPred)){             # if group level predictors at the unit level
       beta.pred <- matrix(beta.pp, ncol=ncol(UnitPred), byrow=TRUE)  # matrix form of the average effects
       UnitCova <- remove.NA(unique(UnitPred))                        # get rid of missing data
       beta.2 <- Random.Coef(Cova=UnitCova, random=random.pred, coef.pred=beta.pred, RBB=RB )
       betaii <- rbind(betaii, beta.2)
     }else{                                    # if no group level predictors
       beta.2 <- matrix(NA, ncol=ncol(random.pred), nrow=GU)  # simply add the random residual on the average
       for (h in 1: GU){
          beta.2[h,] <- random.pred[h,]+beta.pp
       }
      betaii <- rbind(betaii, as.vector(beta.2))
      } 
  } else {                                 # if only random intercept specified, stop grogram and return error message
       warning("No unit-specific effects to return")
       return (NULL)
  }
}

 # Compute time-specific random effects
  if  ("betat"%in% monitor){
       time.random <- as.matrix(new.ct)
      if (length(unique(S))>1){
        if (is.numeric(TimePred)){
          gamma.pp <- new.beta[(ncol(X)-ncol(TimePred)*ncol(S)+1):ncol(X)]
          time.pred <- matrix(gamma.pp, ncol=ncol(TimePred), byrow=TRUE)
          TimeCC <- as.matrix(na.omit((unique(TimePred))))
          gamma.2 <- Random.Coef(Cova=TimeCC, random=time.random, coef.pred=time.pred, RBB=TB)
          gammai <- rbind(gammai, gamma.2)
       }else{
          gamma.pp <- new.beta[(ncol(X)-ncol(S)+1):ncol(X)]
          gamma.2 <- matrix(NA, ncol=ncol(time.random), nrow=GT)
         for (h in 1: GT){
           gamma.2[h,] <- time.random[h,]+gamma.pp
          }
         gammai <- rbind(gammai, as.vector(gamma.2))
       } 
     }else{
        warning("No time-specific effects to return")
        return (NULL)
      }
   }


   #*********************** Update other parameters ************************************
       
       if ("ystar"%in% monitor){
          if ((g%%thin ==0)| g==m){
         ystarm <- rbind(ystarm, new.ystar)
          }
        }
       
         if ("u"%in% monitor){
             if ((g%%thin ==0)| g==m){
               u.aux <- rbind(u.aux, new.u)
             }
         }
       
       if ("bi" %in% monitor){
         bbb <- rbind(bbb, as.vector(new.bi))
       }
       
       if ("ct"%in% monitor){
        timetime <- as.vector(new.ct)
        ccc<-rbind(ccc, timetime)
       }
       
       if ("beta"%in% monitor){
        beta<-rbind(beta,new.beta)
       }
       
       if ("rho"%in% monitor){
         if (lag==1){
           phi <- c(phi,new.phi)
           }else{
           phi <- rbind(phi,new.phi)
          }
       }
       
       if ("D"%in% monitor){
         if (ncol(W)==1){
           DDD <- rbind(DDD, as.vector(Di))
         }else{
           DDD <- rbind(DDD, sD)
         }
       }
       
        if ("E"%in% monitor){
            if (ncol(S)==1){
           EEE <- rbind(EEE, as.vector(Ei))
         }else{
           EEE <- rbind(EEE, sE)
         }
       }
       
       if (categ>2){
        if ("tau"%in% monitor){
        tau <- rbind(tau, new.tau)
       }
      }
       
     }

#__________________________Returning MCMC output_______________________________

  if ("rho"%in% monitor){ 
    if (lag==1){posterior.phi <- phi[-c(1:burnin)]
             }else{
        posterior.phi <- phi[-c(1:burnin),]
     }
  }else{
        posterior.phi <- phi
  }      
  if ("beta"%in% monitor){ 
     posterior.beta <- beta[-c(1:burnin),]
  }else{
     posterior.beta <- beta
   }
  if ("bi"%in% monitor){ 
     posterior.unit.res <- bbb[-c(1:burnin),]
  }else{
     posterior.unit.res <- bbb
  }
 if ("ct"%in% monitor){ 
     posterior.time.res <- ccc[-c(1:burnin),]
  }else{
     posterior.time.res <- ccc
  }   
  if ("E"%in% monitor){ 
     posterior.time.cov <- EEE[-c(1:burnin),]
  }else{
     posterior.time.cov <- EEE
  }    
  if ("D"%in% monitor){ 
     posterior.unit.cov <- DDD[-c(1:burnin),]
  }else{
     posterior.unit.cov <- DDD
  }
  if ("betai"%in% monitor){ 
     posterior.unit <- betaii[-c(1:burnin),]
  }else{
     posterior.unit <- betaii
  }   
  if ("betat"%in% monitor){ 
      posterior.time <- gammai[-c(1:burnin),]
  }else{
      posterior.time <- gammai
  }    
   if ("u"%in% monitor){
      posterior.u.aux <- u.aux[-c(1:(burnin/thin)),]
  }else{
      posterior.u.aux <- u.aux
  }    
   if ("ystar"%in% monitor){ 
     posterior.data <- ystarm[-c(1:(burnin/thin)),]
  }else{
     posterior.data <- ystarm
  }
    if (categ>2){
    if ("tau"%in% monitor){ 
      posterior.tau <- tau[-c(1:burnin),]
     }else{
      posterior.tau <- tau
  }
 }
     if (categ>3){
    return(list(Phi=posterior.phi, Beta=posterior.beta,  bi=posterior.unit.res, ci=posterior.time.res , var.c=posterior.time.cov, Cov.betai=posterior.unit.cov, unit.random=posterior.unit, time.random=posterior.time, cutpoint=posterior.tau, auxiliary.u=posterior.u.aux, data.augmented=posterior.data))
     }else{
    return(list(Phi=posterior.phi, Beta=posterior.beta,  bi=posterior.unit.res, ci=posterior.time.res, var.c=posterior.time.cov, Cov.betai=posterior.unit.cov, unit.random=posterior.unit, time.random=posterior.time, auxiliary.u=posterior.u.aux, data.augmented=posterior.data))

    }
 }


# This code has not been changed completely
#______________ Code for B/OTSCS with Serial Correlation and with U in the same
#_____ Block with beta and {b_i} ____________________________________________________

 CholARP.BTSCSART2 <- function(y, Fixed, UnitRandom, TimeRandom, UnitPred, TimePred, Unit, Time,  m=10000, burnin=5000,  inittau, init.phi, init.ystar=0, init.beta=0,  init.b=0, init.c=0,a1=3, a2=0.1, beta0=rep(0,24),B0=diag(10, 24), D0=diag(1.67, 5) , d0=9, E0=diag(1.67, 2), e0=6,  rescale=1, pgm=1, lag=2){
 
#***************************** Arrange the Data *************************************#
       if(!is.numeric(UnitPred)& !is.numeric(TimePred)){
      NewData <- DataArrange.NP(id1=Unit, id2=Time, y=y, Fixed=Fixed, UnitRandom=UnitRandom, TimeRandom=TimeRandom)
    } else{
        if (!is.numeric(UnitPred) &is.numeric(TimePred)){
           NewData <- DataArrange.TP(id1=Unit, id2=Time, y=y, Fixed=Fixed, UnitRandom=UnitRandom, TimeRandom=TimeRandom, TimePred=TimePred)
        } else{
           if (is.numeric(UnitPred) &!is.numeric(TimePred)){
              NewData <- DataArrange.UP(id1=Unit, id2=Time, y=y, Fixed=Fixed, UnitRandom=UnitRandom, TimeRandom=TimeRandom, UnitPred=UnitPred)
             }else{
               NewData <- DataArrange(id1=Unit, id2=Time, y=y, Fixed=Fixed, UnitRandom=UnitRandom, TimeRandom=TimeRandom, UnitPred=UnitPred, TimePred=TimePred)
             }
        }
   }

  newid1 <- NewData[[1]]
  newid2 <- NewData[[2]]
  Y <- NewData[[3]]
  X <- as.matrix(NewData[[4]])
  W <- as.matrix(NewData[[5]])
  S <- as.matrix(NewData[[6]])
#****************************** Dimensions again ************************************#
   
   # number of all observations
   N <- length(Y)
   # Dimension of all fixed parameter vector beta 
   FB <- ncol(X)
   # Dimension of random fixed parameter vector bi
   RB <- ncol(W)
   # Dimension of random fixed parameter vector ci
   TB <- ncol(S)
 
   # create indices for two non-nested groupings ##########
     entry <- seq(1: N)                   # the global length of observations
     id1 <- Index(newid1)                 # unit id repeating t times t=# observations
     id2 <- (newid2)-min(newid2)+1        # year id repeating n times n=# unit in each time period
     GU <- max(id1)                       # how many units
     GT <- max(id2)                       # how many time periods in total
     categ <- length(unique(Y))           # binary or ordinal
   
   # output the summary of observations
      sample.size.unit <- as.vector (table (id1))
      sample.size.time <- as.vector (table (id2))
      cat("number of observations for each unit: ",  sample.size.unit,"\n", spe=" ")
      cat("number of observations for each time period: ",  sample.size.time,"\n", sep=" ")
 
#**************************** initial values **************************************#
   
     if (lag==1){
        phi <- init.phi
      }else{
        phi <- matrix(init.phi, nrow=1, ncol=lag)
         }
      beta <- matrix(init.beta, nrow=1, ncol=FB)
      bbb <- matrix(init.b, nrow=1, ncol=RB*GU)
      DDD <- matrix(0, nrow=1, ncol=(RB^2-RB)/2+RB)
      ccc <- matrix(init.c, nrow=1, ncol=GT*TB)
      EEE <- matrix(0, nrow=1, ncol=(TB^2-TB)/2+TB)
      betaii <- matrix(0, nrow=1, ncol=RB*GU)
      gammai <- matrix (0, nrow=1, ncol=GT*TB)
      ystarm <- matrix(init.ystar, nrow=1, ncol=N)
      u.aux <-  matrix(rnorm(N), nrow=1, ncol=N)
       if (categ>2){
        tau <- matrix(inittau, nrow=1, ncol=(categ-2))
          }
   
############### loop begins ################################################

 #Tracking the process of MCMC simulation
    for(g in 2:(m+burnin)) { 
       if ((g%%5 ==0)| g==m){
          cat("g=",g,"\n")
        }
       
# Define variables which will be used throughout the iteration 
 if (g==2){
   current.bi <- matrix(init.b, ncol=RB)           # current values of bi
   current.ct <- matrix(init.c, ncol=TB)           # current values of ct
   current.beta <- init.beta                         # current values of beta
   current.ystar <- init.ystar                     # the current value of ystar
   current.phi <- init.phi
   current.u <- init.u
   if (categ>2){
     current.tau <- init.tau
   }
 }else{
   current.bi <- new.bi                            # current values of bi
   current.ct <- new.ct                            # current values of ct
   current.beta <- new.beta                        # current values of beta
   current.ystar <- new.ystar                      # the current value of ystar
   current.phi <- new.phi
   current.u <- new.u
  if (categ>2){
     current.tau <- new.tau
   }
 }
  

 #_______________Block 1: variance Matrix of E _____________________________

        OEin <- solve(E0)+sumtrans(Mat=current.ct)
        sE<-rWishart (1, e0+GT, solve(OEin))
        Ei <- Wishart(element=sE, col=TB)
       
#_______________Block 2: variance Matrix of D _____________________________
       
       ODin <- solve(D0)+sumtrans(Mat=current.bi)
       sD<-rWishart (1, d0+GU, solve(ODin))
       Di <- Wishart(element=sD, col=RB)
 #__________________________________________________________________________

   if (lag==1){
      var.phig <- 1/(1- current.phi^2)                             # the variance computed by using the current values of phi
      TCov <- Cov.AR1(col=GT, rho= current.phi,sigma=var.phig)
   }else{
       CCov <- Cov.ARp(phig=current.phi, col=GT)        
       TCov <- CCov[[2]]
       Sig <- TCov[1:lag, 1:lag]
       SSig<-solve(Sig) 
       var <- Sig[1,1]
      }
       
##### Variables for outputs from the loop ################################
       
      betaterm1 <- matrix(0, nrow=FB, ncol=FB)    # for posterior covariance of beta
      betaterm2 <- rep(0, FB)                     # for posterior mean of beta
 
      new.ystar <- rep(NA, N)

   for (n in 1: GU){                                # loop for each unit
      rows <- entry[which(id1==n)]                    # the obserbations of unit n
      nT <- length(rows)
      yi <- Y[rows]
      ystari <- current.ystar[rows]
      Xi <- X[rows,]
      Wi <- W[rows,]
      Si <- S[rows,]
      uuu <- current.u[rows]
      bbbn <- current.bi[n,]
      TT <- id2[rows]
      Tccin <- current.ct[TT,]
     
       rant <- c()
         for (l in 1:nT){
          rant[l] <- Si[l,]%*%cbind(Tccin[l,])
        }

   #__________________________Block 3:Updata U _____________________________
     
        nCov <- TCov[1:nT, 1:nT]
        nVK <- V.matrix(Omega=nCov, col=nT)
        nV <- nVK[[2]]                          
        nK <- nVK[[1]]
        u.decom <- nV%*%cbind(uuu)
 
  #__________________________Block 4: Updata ystar ________________________
     
         mean.ystar <- Xi%*% cbind(current.beta)+Wi%*%cbind(bbbn)+rant+u.decom
         ystar <- c()
          if (categ==2){
            for (r in 1: nT){
               ystar[r] <- BTnormlog(x=yi[r], mu=mean.ystar[r], sigma=sqrt(nK))
               }
          }else{
            for (j in 1:nT){
               ystar[j] <- OrdinalAugment(obs=yi[j], mu=mean.ystar[j], sigma=sqrt(nK), maxcat=max(Y), mincat=min(Y), cutpoint=current.tau)
           }
         }
        new.ystar[rows] <- ystar
        
   #_________________ calculate BB1 and BB2 for beta __________________________
   
      var.beta <- pseudoinverse(nCov+Wi%*%solve(Di)%*%t(Wi))
      betaterm1 <- betaterm1+t(Xi) %*% var.beta %*%Xi
      betaterm2 <- betaterm2+t(Xi)%*% var.beta %*% (ystar-rant)
         }
       
 # End of the loop
 
 #_______________Block 5 Beta and {b_i}: ______________________________________
       # 1) \beta
       B1<-pseudoinverse(betaterm1+solve(B0))
       beta1<-t(B1%*%(betaterm2+solve(B0)%*%beta0))
       #B1 <- make.positive.definite(B1)
       new.beta<-mvrnorm(1, beta1, B1)
    #   cat("beta ", sum(new.beta), "\n", sep=" ")

       # 2) {b_i}
      new.bi <- matrix(NA, nrow=1, ncol=RB)
      new.u<- rep(NA, N)
      kap <- rep(NA, GU)                                 # Store the kappa for variance
      qu <- rep(NA,N)                             # Store the output of the auxiliary term

       for (n in 1: GU){                                # loop for each unit
         rows <- entry[which(id1==n)]                    # the obserbations of unit n
         nT <- length(rows)
         nCov <- TCov[1:nT, 1:nT]
         nVK <- V.matrix(Omega=nCov, col=nT)
         nV <- nVK[[2]]                          
         nK <- nVK[[1]]
         ystari <- new.ystar[rows]
         Xi <- X[rows,]
         Wi <- W[rows,]
         Si <- S[rows,]
         TT <- id2[rows]
         Tccin <- current.ct[TT,]
         nkap <- kap[n]
     
         rant <- c()
         for (l in 1:nT){
          rant[l] <- Si[l,]%*%cbind(Tccin[l,])
          }
     
      Cov.b <- pseudoinverse(Di + t(Wi)%*%solve(nCov)%*%Wi)
      mean.b <- Cov.b%*% t(Wi)%*%solve(nCov)%*%(ystari-Xi%*% cbind(new.beta)-rant)
      #Cov.b <- make.positive.definite(Cov.b)
      b.drawn <- mvrnorm(1, mu=mean.b, Sigma=Cov.b)
      new.bi <- rbind(new.bi, b.drawn)
         
        Cov.U <- solve(diag(1, nT) + 1/nK*t(nV)%*%nV)   
        mean.U <- Cov.U %*% t(nV) %*% (ystari-Xi%*%cbind(new.beta)-Wi%*%cbind(b.drawn)-rant)/nK
        u.draw <- mvrnorm(1, mu=mean.U, Sigma=Cov.U)
        new.u[rows] <- u.draw
        Nu.decom <- nV%*%cbind(u.draw)
        qu[rows] <- Nu.decom
        kap[n] <- nK
   }

        new.bi <- new.bi[-1,]
       
  #___________________ Block 6: Time-Specific Random Effect ________________

         new.ct <-matrix(NA, nrow=1, ncol=TB)
        for (j in 1: GT){
            row1<-entry[which(id2==j)]  # Those entries with t=j
            H <- length(row1)           #  #observations at t period
            iden1<-id1[row1]            # The units at t period
            tttn <- new.bi[iden1,]         # The random unit-specific effects of
            yt1<-new.ystar[row1]
            Xt1<-X[row1,]
            Wt1<-W[row1,]
            St1 <- S[row1,]
            Ut1 <- qu[row1]
            kapt1 <- kap[iden1]
            
            term7 <- c()
          for (q in 1:H){
            term7[q] <- Wt1[q,]%*% cbind(tttn[q,])
              }
            
             PM <- matrix(0, nrow=H, ncol=H)
          for (c in 1:H){
             PM[c,c] <- kapt1[c]
              }
            
           cov.S <- solve(Ei+t(St1)%*%solve(PM)%*%St1)
           mean.S <- cov.S%*% t(St1) %*%solve(PM)%*%(yt1-term7-Xt1%*%cbind(new.beta)-Ut1)
           supd <- rmvnorm(1, mean.S, cov.S)
          
           new.ct <- rbind(new.ct, supd)
             }
       
           new.ct <- new.ct[-1,]
#__________________Block 7: Autoregressive Coef.s_________________________

     if (lag==1){
       PhiPr <- PhiTerm0(Entry=entry, unit=id1, NN=GU, yystar=new.ystar,newb=new.bi, PX=X, PW=W, PS=S, Ti=id2, CC=new.ct, newbeta=new.beta)
       PhiSig <- PhiPr[1]
       PhiMean <- PhiPr[2]
       hatPhi<-1/PhiSig
       hatphi<-hatPhi*PhiMean
       phi.draw <- rtnorm(1, hatphi, sqrt(hatPhi), lower=-1, upper=1)
  cat("phi.draw ", phi.draw, "\n", sep=" ")
        ####### Accept and Reject ########
       
       alpha<-AcceptRate0(rho=phi.draw, orho=phig, unit=id1, precision=1,AY=new.ystar, AX=X, AW=W, AS=S, Abeta=new.beta, Ab=new.bi, AT=id2, Ac=new.ct, Entry=entry)

        u <- runif(1)
        if( u <= alpha ) {
         new.phi<-phi.draw
           }else{
            phi<-current.phi
           }
       }
      if (lag>1){ 
         PhiPr <- PhiTermARpT(Entry=entry, unit=id1, NN=GU, yystar=new.ystar, newb=new.bi, PX=X, PW=W, PS=S, Ti=id2, CC=new.ct, newbeta=new.beta, lag=lag)
         PhiSig <- PhiPr[,-(lag+1)]
         PhiMean <- PhiPr[,(lag+1)]
         phigg<-PhiPropARp(Ymt=PhiSig, Ymt2=PhiMean,scale=1, Oldphi=phig)
         cat("phi.draw ", phigg, "\n", sep=" ")
         
 ####### Accept and Reject ########
         
        alpha<-AcceptRateARpT(nrho=phigg, orho=phig, Osig=Sig , unit=id1, AY=new.ystar, AX=X, AW=W, AS=S, Abeta=new.beta, Ab=new.bi,  AT=id2, Ac=new.ct, Entry=entry, f=g, lag=lag)
     #   cat("alpha ", alpha, "\n", sep=" ")
       u <- runif(1)
        if( u <= alpha ) {
            new.phi <- phi.draw
           } else {
            new.phi<- current.phi
             }
          }
  #________________Block8 _Cutpoints Updating ________________________#

       if (categ>2){
       new.tau <- CutPointUpdate(obs=Y, latent=new.ystar, oldcut=current.tau)
          }

 # __________________________ PGM Updating_________________________________
       
       if (pgm==1){
           if (any(beta0!=0)){
        warning("priors have to be centered at 0 for using PGM-MGMC updating")
        return(NULL)
         }
             phig <-new.phi   
             if (lag==1){
             Phig <- 1/(1-phig^2)
              TCov <- Cov.AR1(col=GT, rho=phig,sigma=Phig)
                }else{
               CCov <- Cov.ARp(phig=phig, col=GT)        
               TCov <- CCov[[2]]
               }
             qq1 <- 0
             qq2 <- 0
          for (f in 1: GU){
            row2 <- entry[which(id1==f)]    # Those entries with n=i
            H <- length(row2)               
            iden3 <-id2[row2]               # The units at t period
            tttn <- new.bi[f,]                 # The random unit-specific effects
            yt1 <- new.ystar[row2]
            Xt1 <- X[row2,]
            Wt1 <- W[row2,]
            St1 <- S[row2,]
            nCov <- TCov[1:H, 1:H]
            timec <- new.ct[iden3,]
               rrant <- c()
              for (w in 1:H){
                rrant[w] <- St1[w,]%*% cbind(timec[w,])
                }
             termm <- yt1-Xt1%*%as.matrix(new.beta)-Wt1%*%as.matrix(tttn)-rrant
             qq1 <- qq1+t(termm)%*%solve(nCov)%*%termm
             qq2 <- qq2+rbind(tttn)%*%Di%*%cbind(tttn)
            }
           
             qq3 <- rbind(new.beta)%*% solve(B0)%*% cbind(new.beta)
           
             qq4 <- 0
             for (t in 1: GT){
               qq4 <- qq4+rbind(new.ct[t,])%*%Ei%*%cbind(new.ct[t,])
               }
           
             bb <- (qq1+qq2+qq3+qq4)/2
                         
           if (categ>2){
               aa <- (N+ ncol(beta)+ncol(bbb)+ncol(ccc)+length(new.tau)+1)/2
                  } else{
               aa <- (N+ ncol(beta)+ncol(bbb)+ncol(ccc)+1)/2
           }

           gamm <- sqrt(rgamma(1, shape=aa, scale=1/bb))
           new.ystar <- new.ystar*gamm
           new.beta <- new.beta*gamm
           new.bi <- new.bi*gamm
           new.ct <- new.ct*gamm
           
           if (categ>2){
           new.tau <- new.tau*gamm
             }
            cat("gamm ", gamm, "\n", sep=" ")
         }
   #_____________________Random Effects ____________________________________
  if ("betai"%in% monitor){
# Unit-specific random effects
   random.pred <- as.matrix(new.bi)           # the new draw of random unit effects
   beta.pp <- beta[g,][(ncol(Fixed)+1):(ncol(X)-TB*ncol(S))]
   beta.pred <- matrix(beta.pp, ncol=ncol(UnitPred), byrow=TRUE)
   UnitCova <- remove.NA(unique(UnitPred))
   beta.2 <- Random.Coef(Cova=UnitCova, random=random.pred, coef.pred=beta.pred, RBB=RB )
   betaii <- rbind(betaii, beta.2)
 }
  if  ("betat"%in% monitor){
# Time-specific random effects
   time.random <- as.matrix(new.ct)
   gamma.pp <- beta[g,][(ncol(X)-TB*ncol(S)+1):ncol(X)]
   time.pred <- matrix(gamma.pp, ncol=ncol(TimePred), byrow=TRUE)
   TimeCC <- as.matrix((unique(TimePred))[-1,] )  
     gamma.2 <- Random.Coef(Cova=TimeCC, random=time.random, coef.pred=time.pred, RBB=TB)
     gammai <- rbind(gammai, gamma.2)
 }

      if ("ystar"%in% monitor){
          if ((g%%thin ==0)| g==m){
         ystarm <- rbind(ystarm, new.ystar)
          }
        }
         if ("u.aux"%in% monitor){
             if ((g%%thin ==0)| g==m){
               u.aux <- rbind(u.aux, new.u)
          }
       }
       if ("bi" %in% monitor){
         bbb <- rbind(bbb, as.vector(new.bi))
       }
       if ("ct"%in% monitor){
        timetime <- as.vector(new.ct)
        ccc<-rbind(ccc, timetime)
       }
       if ("beta"%in% monitor){
        beta<-rbind(beta,new.beta)
       }
       if ("rho"%in% monitor){
         if (lag==1){
           phi <- c(phi,new.phi)
           }else{
           phi <- rbind(phi,new.phi)
          }
       }
       if ("D"%in% monitor){
        DDD <- rbind(DDD, sD)
       }
        if ("E"%in% monitor){
        EE <- rbind(EEE, sE)
       }
       if (categ>2){
        if ("tau"%in% monitor){
        tau <- rbind(tau, new.tau)
       }
      }
       
     }

     
#__________________________Returning MCMC output_______________________________
   
    if (lag==1){posterior.phi <- phi[-c(1:burnin)]
             }else{
        posterior.phi <- phi[-c(1:burnin),]
     }
   
  posterior.beta <- beta[-c(1:burnin),]
  posterior.unit.res <- bbb[-c(1:burnin),]
  posterior.time.res <- ccc[-c(1:burnin),]
  posterior.time.cov <- EEE[-c(1:burnin),]
  posterior.unit.cov <- DDD[-c(1:burnin),]
  posterior.unit <- betaii[-c(1:(burnin/thin)),]
  posterior.time <- gammai[-c(1:(burnin/thin)),]
  posterior.u.aux <- u.aux
  posterior.data <- ystarm
   if (categ>3){posterior.tau <- tau[-c(1:burnin),]}
       
     if (categ>3){
    return(list(Phi=posterior.phi, Beta=posterior.beta,  bi=posterior.unit.res, ci=posterior.time.res , var.c=posterior.time.cov, Cov.betai=posterior.unit.cov, unit.random=posterior.unit, time.random=posterior.time, cutpoint=posterior.tau, auxiliary.u=posterior.u.aux, data.augmented=posterior.data))
     }else{
    return(list(Phi=posterior.phi, Beta=posterior.beta,  bi=posterior.unit.res, ct=posterior.time.res, var.c=posterior.time.cov, Cov.betai=posterior.unit.cov, unit.random=posterior.unit, time.random=posterior.time, auxiliary.u=posterior.u.aux, data.augmented=posterior.data))

    }
 }
